\[\textbf{Lineu Alberto Cavazani de Freitas}\] \[\textbf{Prof. Cesar Augusto Taconeli}\] \[\textbf{Prof. José Luiz Padilha da Silva}\]

Análise Comportamental de Ovelhas Submetidas a Intervenção Humana usando GAMLSS

Material Rbras: Número de Mudanças de Postura de Orelha e Cauda

1 Experimento

O conjunto de dados é proveniente de um experimento sobre o comportamento de ovelhas, conduzido na fazenda experimental INRA La Fage, Roqueford, França, em setembro de 2015 com o objetivo de verificar o efeito de linhagem genética, escovação e isolamento nas respostas comportamentais dos animais (TAMIOSO et al., 2017). Na ocasião, vinte ovelhas classificadas como reativas ou não reativas ao isolamento social temporário foram submetidas à escovação por um humano familiar. As ovelhas tinham 15 meses de idade, eram não gestantes e não amamentavam quando foram observadas.

O experimento foi conduzido em três sessões experimentais: na primeira tinha-se uma grade de metal separando o animal testado dos demais animais, sem distância entre eles. Na segunda havia duas grades de metal separando os animais a uma distância de 1,7 metros, ou seja, foi imposta a condição de isolamento social. E na terceira sessão, os animais voltaram a ser separados por apenas uma grade:

As sessões de testes ocorreram dois dias após a fase de adaptação dos animais ao equipamento e aos humanos e, em cada sessão, as ovelhas foram observadas em 3 momentos distintos: fase de pré escovação, com duração de 2 minutos e 30 segundos; fase de escovação, com duração de 3 minutos; e pós escovação, com duração de 2 minutos e 30 segundos. O delineamento para um animal pode ser represntado da seguinte forma:

Os dados coletados dizem respeito ao número de mudanças de postura dos animais e à proporção do tempo em que os animais permaneceram em determinadas posturas, tratando-se então, de um conjunto de dados com múltiplas respostas em que não há observações independentes, já que cada animal contribui com nove medidas. Portanto, há a necessidade de incorporar as correlações entre as medidas num mesmo animal e do animal dentro de cada sessão experimental, além da correlação entre as respostas.

Foram avaliados os efeitos de:

Sessão: Fator de 3 níveis que indica a sessão experimental (Se1, Se2, Se3).
Momento: Fator de 3 níveis em que indica o momento experimental (antes, durante, depois).
Linhagem: Fator de 2 níveis que classifica os animais como reativos ou não reativos ao isolamento social temporário.

2 Análise Exploratória

2.1 Dados

setwd("~/Área de Trabalho/PESQUISA/IC_nova/1. Dados")

## Pacotes
library(ggplot2)
library(gridExtra)
library(gamlss)
library(knitr)

######################################################################
## Dados

count <- read.csv2("DataCont.csv", header = T, sep = ',')

## Alterando os contrastes default para os contrastes de interesse

count$sessaoR <- count$sessao

count$sessaoR <- relevel(count$sessaoR, ref = 'Se3')
contrasts(count$sessaoR) <- 'contr.helmert'

contrasts(count$sessaoR)[,1] <- contrasts(count$sessaoR)[,1]/2
contrasts(count$sessaoR)[,2] <- contrasts(count$sessaoR)[,2]/3

## Criando um fator combinando animal e sessao
count$animals <- factor(count$animal):count$sessao

2.2 Histograma e Box-plots

2.3 Gráficos de Perfil

3 Modelos para o Número de Mudanças de Postura de Orelha

Para o ajuste dos modelos de regressão foram consideradas as distribuições:

PO: Poisson
NBI:Binomial Negativa Tipo I
ZIP: Poisson Inflacionada de Zeros
ZINBI: Binomial Negativa Inflacionada de Zeros
ZAP: Poisson Zero Ajustada
ZANBI: Binomial Negativa Zero Ajustada

Para cada uma das 7 distribuições ajustou-se o modelo no qual era considerado sessão, momento, linhagem, as interações duas a duas e efeitos aleatórios a nível de animal e animal dentro de sessão para a média.

Nas distribuições com parâmetros referentes à inflação em 0 também foram acrescentadas as mesmas covariáveis, com exceção dos efeitos aleatórios.

Todas análises foram feitas utilizando o pacote gamlss do software estatístico R.

modelo_norelha <- function(familia){
  gamlss(norelha ~ sessaoR + tempo + linhagem 
         + (sessaoR + tempo + linhagem)^2 +
           re(random = list(animal=~1, animals=~1), method = 'REML'),
         nu.formula = norelha ~ sessaoR + tempo + 
           linhagem+ (sessaoR + tempo + linhagem)^2,
         data = count, method = RS(50), 
         family = familia)}

mPO    <- modelo_norelha(familia = PO())
mNBI   <- modelo_norelha(familia = NBI())
mZIP   <- modelo_norelha(familia = ZIP())
mZINBI <- modelo_norelha(familia = ZINBI())
mZAP   <- modelo_norelha(familia = ZAP())
mZANBI <- modelo_norelha(familia = ZANBI())

3.1 Seleção da Família Via Medidas de Qualidade

A escolha da distribuição utilizada baseou-se nos valores observados para as medidas AIC, BIC, Deviance, Verossimilhança, graus de liberdade e por fim, verificou-se se houve convergência.

Modelo	AIC	BIC	Deviance	logLik	df	Conv	Iter
PO	1191.028	1344.848	1094.6787	-547.3393	48.17479	TRUE	5
NBI	1054.742	1206.695	959.5615	-479.7808	47.59021	TRUE	7
ZIP	1127.398	1315.057	1009.8519	-504.9260	58.77281	TRUE	17
ZINBI	2488.164	2687.903	2363.0517	-1181.5259	62.55619	FALSE	50
ZAP	1137.740	1324.007	1021.0663	-510.5332	58.33675	TRUE	12
ZANBI	1046.170	1247.390	920.1308	-460.0654	63.01976	TRUE	10

A distribuição selecionada foi a Binomial Negativa Zero Ajustada.

3.2 Reajustes

Considerando o modelo ZANBI, foram realizados 3 reajustes para o modelo original:

Retirando interações de \(\nu\)
Retirando interações de \(\mu\) e \(\nu\)
Retirando interações de \(\mu\) e todas as covariáveis de \(\nu\)

Todos os modelos reajustados são modelos encaixados ao modelo original, sendo assim utilizou-se o teste de razão de verossimilhanças para selecionar o menor modelo possível que não difere estatísticamente do original.

3.2.1 Retirando interações de \(\nu\)

mZANBIr1 <- gamlss(formula = norelha ~ sessaoR + tempo + linhagem + 
                  (sessaoR + tempo + linhagem)^2 +
                  re(random = list(animal = ~1, animals = ~1), method = 'REML'), 
                nu.formula = norelha ~ sessaoR + tempo + linhagem, 
                family = ZANBI(), data = count)

3.2.2 Retirando interações de \(\mu\) e \(\nu\)

mZANBIr2 <- gamlss(formula = norelha ~ sessaoR + tempo + linhagem + 
                  re(random = list(animal = ~1, animals = ~1), method = 'REML'), 
                nu.formula = norelha ~ sessaoR + tempo + linhagem, 
                family = ZANBI(), data = count)

3.2.3 Retirando interações de \(\mu\) e todas as covariáveis de \(\nu\)

mZANBIr3 <- gamlss(formula = norelha ~ sessaoR + tempo + linhagem + 
                  re(random = list(animal = ~1, animals = ~1), method = 'REML'),
                family = ZANBI(), data = count)

3.2.4 Medidas de qualidade e TRV

##       Modelo      AIC      BIC Deviance    logLik       df Conv P-val
## 1   Completo 1046.170 1247.390 920.1308 -460.0654 63.01976 TRUE     -
## 2 Reajuste 1 1036.359 1212.034 926.3192 -463.1596 55.01976 TRUE  0.63
## 3 Reajuste 2 1022.671 1174.258 927.7194 -463.8597 47.47563 TRUE  0.95
## 4 Reajuste 3 1037.248 1172.871 952.2970 -476.1485 42.47563 TRUE  0.05

Conclusões

O modelo sem interações para nu não difere do modelo completo.
O modelo sem interações para mu e nu não difere do modelo completo.
O modelo sem interações para mu e nenhuma covariável para nu difere do modelo com covariáveis para nu.
Sendo assim, as covariáveis para mu e nu devem ser mantidas, porém sem as interações.

4 Modelos para o Número de Mudanças de Postura de Cauda

Para o número de mudanças de postura de cauda foi experimentado um modelo com os efeitos fixos de sessão, momento e linhagem na média e na dispersão. Além disso considerou-se o efeito aleatório a nível de animal na média.

mNBI2_sigma <- 
  gamlss(ncauda ~ sessaoR + tempo + linhagem +
           re(random = list(animal=~1)),
         sigma.formula = ~  sessaoR +tempo+ linhagem,
         data = count, family = NBII(), method = RS(100))

5 Análise de Resíduos

Parte importante da etapa de ajuste de modelos de regressão é a análise de resíduos. A análise de resíduos em modelos da classe GAMLSS é baseada em resíduos quantílicos aleatorizados, estes resíduos são de interpretação mais simples visto que, sob um modelo bem ajustado, possuem distribuição Normal.

6 Resumo dos modelos finais

6.1 Número de Mudanças de Postura de Orelha

summary(mZANBIr2)

## ******************************************************************
## Family:  c("ZANBI", "Zero Altered Negative binomial type I") 
## 
## Call:  gamlss(formula = norelha ~ sessaoR + tempo + linhagem +  
##     re(random = list(animal = ~1, animals = ~1), method = "REML"),  
##     nu.formula = norelha ~ sessaoR + tempo + linhagem,  
##     family = ZANBI(), data = count) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.49220    0.07710  32.324  < 2e-16 ***
## sessaoR1      0.42218    0.10556   3.999 0.000105 ***
## sessaoR2     -0.16399    0.09251  -1.773 0.078580 .  
## tempoDepois  -0.43246    0.09736  -4.442 1.86e-05 ***
## tempoDurante -1.05623    0.11387  -9.276 4.36e-16 ***
## linhagems+    0.28393    0.08580   3.309 0.001205 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -1.9154     0.2149  -8.912 3.43e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  logit 
## Nu Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -4.1780     1.0714  -3.899 0.000153 ***
## sessaoR1      -2.6613     1.0911  -2.439 0.016049 *  
## sessaoR2       0.4330     0.7312   0.592 0.554710    
## tempoDepois    1.4892     1.1487   1.296 0.197069    
## tempoDurante   2.7598     1.0831   2.548 0.011973 *  
## linhagems+    -1.1191     0.6379  -1.754 0.081693 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms maybe are not accurate. 
## ------------------------------------------------------------------
## No. of observations in the fit:  180 
## Degrees of Freedom for the fit:  47.47563
##       Residual Deg. of Freedom:  132.5244 
##                       at cycle:  14 
##  
## Global Deviance:     927.7194 
##             AIC:     1022.671 
##             SBC:     1174.258 
## ******************************************************************

6.2 Número de Mudanças de Postura de Cauda

summary(mNBI2_sigma)

## ******************************************************************
## Family:  c("NBII", "Negative Binomial type II") 
## 
## Call:  gamlss(formula = ncauda ~ sessaoR + tempo + linhagem +  
##     re(random = list(animal = ~1)), sigma.formula = ~sessaoR +  
##     tempo + linhagem, family = NBII(), data = count,  
##     method = RS(100)) 
## 
## Fitting method: RS(100) 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.03866    0.14764   7.035 6.35e-11 ***
## sessaoR1      0.07220    0.17355   0.416    0.678    
## sessaoR2      0.05521    0.14472   0.382    0.703    
## tempoDepois  -0.09011    0.21369  -0.422    0.674    
## tempoDurante  1.66707    0.15827  10.533  < 2e-16 ***
## linhagems+    0.16417    0.13724   1.196    0.233    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   0.80294    0.38520   2.084   0.0388 * 
## sessaoR1     -0.56679    0.42998  -1.318   0.1894   
## sessaoR2     -0.07557    0.37539  -0.201   0.8407   
## tempoDepois   0.55825    0.49029   1.139   0.2567   
## tempoDurante  1.15304    0.42972   2.683   0.0081 **
## linhagems+   -0.00657    0.36109  -0.018   0.9855   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms maybe are not accurate. 
## ------------------------------------------------------------------
## No. of observations in the fit:  180 
## Degrees of Freedom for the fit:  27.95139
##       Residual Deg. of Freedom:  152.0486 
##                       at cycle:  13 
##  
## Global Deviance:     913.2412 
##             AIC:     969.1439 
##             SBC:     1058.392 
## ******************************************************************

TAMIOSO, P. R.; BOISSY, A.; BOIVIN, X.; CHANDEZE, H.; ANDANSON, S.; DELVAL, E.; HAZARD, D.; TACONELI, C. A.; MOLENTO, C. F. M. Does emotional reactivity influence behavioral and cardiac responses of ewes submitted to brushing? Behavioural Processes, p. np, 2017.